program main_st_lrtest_bootstrap
use studentscr
use matrixfuns
use datahadler
use stestmod2
use ghlliks
implicit none
integer, parameter:: TT=1000,nn=3,S=10000,ngrid=10,sboot=100
integer,parameter:: nvar0=2*nn+6
integer, parameter:: nvar=nvar0+nn+1
integer iparam(7),T,n,i,j,sit,seedvec(S),seedind(2),seedboot(sboot),inn,id

double precision:: y(TT,nn),ytot(TT,nn),yi(TT,nn)&
				 &,thscale(nvar),thiterfix(nvar),rparam(7),a1,a2,c9

double precision:: thlb(nvar),thub(nvar),thtrue(nvar),th0(nvar)&
				 &,th(nvar),llik,grad(nvar),therr(nvar),lika,likb,tha(nvar),thb(nvar)

double precision:: mth(S,nvar),lr(S)

double precision:: sigrid(ngrid)&
&,llikf,llikj,thfix(nvar),llik1,matthshape(S,nn+2)&
&,lri(sboot),pvali(sboot),plrboot(S)

double precision:: eta_0,grd1(nvar),grd2(nvar)

double precision, parameter:: pi=3.14159265358979d0
				 
common T,n,y
common /eta_prev_step/ eta_0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
T=TT
n=nn
call rnset(1)
!!!!!!!!!!True parameters
a1=.1d0
a2=.85d0
c9=.999d0

thtrue(1:n)=vassign(n,.2d0)			!delta
thtrue(n+1:2*n-1)=vassign(n-1,1.0d0)	!c
thtrue(2*n)=0.05d0					!alpha0
thtrue(2*n+1)=0.05d0				!phi0
thtrue(2*n+2)=dasin(dsqrt(a1/c9))    !alpha1
thtrue(2*n+3)=dasin(dsqrt(a2/(c9-a1)))  !alpha2
thtrue(2*n+4)=dasin(dsqrt(a1/c9))		!phi1
thtrue(2*n+5)=dasin(dsqrt(a2/(c9-a1)))	!phi2
thtrue(2*n+6)=.1d0
!!!!!!!!! Limits

thlb(1:n)=vassign(n,-1.0d1)			   !delta
thlb(n+1:2*n-1)=vassign(n-1,-1.0d1)		   !c
thlb(2*n)=1.0d-3					   !alpha0
thlb(2*n+1)=1.0d-3					   !phi0
thlb(2*n+2:2*n+5)=vassign(4,0.0d0)
thlb(2*n+6)=1.0d-3
thlb(2*n+7)=1.0d-1
thlb(2*n+8:3*n+7)=vassign(n,-100.0d0)



thub(1:n)=vassign(n,1.0d1)			   !delta
thub(n+1:2*n-1)=vassign(n-1,1.0d1)		   !c
thub(2*n)=10.0d0 					   !alpha0
thub(2*n+1)=10.0d0					   !phi0
thub(2*n+2:2*n+5)=vassign(4,2.0d0*pi)
thub(2*n+6)=0.499d0
thub(2*n+7)=0.99d0
thub(2*n+8:3*n+7)=vassign(n,100.0d0)


thlb(2*n+2:2*n+5)=vassign(4,1.0d-2)
thub(2*n+2:2*n+5)=vassign(4,2.0d0*pi-1.0d-2)
thlb(2*n+3)=.1d0
thlb(2*n+5)=.1d0


thscale=vassign(nvar,1.0d0)

sigrid=linspace(8.0d-1,.99d0,ngrid)


call readvecint(seedvec,S,"seedsimst.txt",1)
call readvecint(seedind,2,"seedind.txt",1)


do sit=seedind(1),seedind(2),1
print*, sit

call rnset(seedvec(sit))

!!!!!!!!!!Simulation
call stsimul(nvar,T,n,thtrue,ytot)


y=ytot(1:T,1:n)
!goto 10




call erset(0,0,0)

!!!!!!!!!!Estimation
call du4inf(iparam,rparam)
iparam(3)=1000
iparam(4)=1000
iparam(5)=1000
!rparam(1)=1.0d-7
!rparam(2)=1.0d-15
call dbcong(stllik,stscr,nvar,thtrue,0,thlb,thub&
          &,thscale,1.0d0,iparam,rparam,th,llik)

thiterfix=th
therr=dabs(thub-th)*dabs(th-thlb)
call dsvrgn(nvar0,therr(1:nvar0),therr(1:nvar0))

mth(sit,:)=th
if (therr(1)>0.0d0) then
	th(nvar0+1)=0.99d0
	th(nvar0+2:nvar)=vassign(n,0.1d0)
	th0=th

	call stllik(nvar0,th(1:nvar0),llikf)

	do i=1,ngrid,1
		th(nvar0+1)=sigrid(i)
		call ghllik(nvar,th,llikj)
		if (llikj<llikf) then
			th0=th
			llikf=llikj
		end if
	end do
	!!GRID
	call erset(0,0,0)
	call du4inf(iparam,rparam)
	iparam(3)=100
	iparam(4)=1000
	iparam(5)=1000

	call dbcong(ghllik,ghscr,nvar,th0,0,thlb,thub&
          &,thscale,1.0d0,iparam,rparam,th,llik1)

	lr(sit)=-2.0d0*(llik1-llik)
    matthshape(sit,:)=th(nvar0:nvar)
	print*,sit,lr(sit)
	call exportvec(th,"thint.txt",15)
	
	do j=1,sboot,1
	print*,j
!	call rnget(id)
!	seedboot(j)=id
!	call exportvecint(seedboot,"datseed.txt",20)
!print*, "seed",id
		call stsimul(nvar,T,n,thiterfix(1:nvar0),yi)

		y=yi(1:T,1:n)

		!!!!!!!!!!!!!!
		call du4inf(iparam,rparam)
		iparam(3)=100
		iparam(4)=1000
		iparam(5)=1000
		call dbcong(stllik,stscr,nvar,thiterfix,0,thlb,thub&
          &,thscale,1.0d0,iparam,rparam,th,llik)

		th(nvar0+1)=0.99d0
		th(nvar0+2:nvar)=vassign(n,0.01d0)
    	th0=th

		call stllik(nvar0,th(1:nvar0),llikf)

		do i=1,ngrid,1
			th(nvar0+1)=sigrid(i)
			call ghllik(nvar,th,llikj)
			if (llikj<llikf) then
				th0=th
				llikf=llikj
			end if
		end do
		!!GRID
		call erset(0,0,0)
		call du4inf(iparam,rparam)
		iparam(3)=100
		iparam(4)=1000
		iparam(5)=1000
		call dbcong(ghllik,ghscr,nvar,th0,0,thlb,thub&
          &,thscale,1.0d0,iparam,rparam,th,llik1)

		lri(j)=-2.0d0*(llik1-llik)
	end do

	where (lri .gt. lr(sit))
		pvali=1.0d0
	elsewhere
		pvali=0.0d0
	end where
		plrboot(sit)=sum(pvali)/dble(sboot)
end if

	call exportvec(vec2(mth(seedind(1):sit,:)),"matth2.txt",1)
	call exportvec(vec2(matthshape(seedind(1):sit,:)),"matthsb.txt",1)
	call exportvec(lr(seedind(1):sit),"lrtests.txt",1)
	call exportvec(plrboot(seedind(1):sit),"pvlrstb.txt",1)

end do

end program main_st_lrtest_bootstrap